In [1]:
from __future__ import division
import numpy as np
import pandas as pd
from scipy import stats
import statsmodels.formula.api as smf
import statsmodels.stats.api as sms
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
Давайте проанализируем данные опроса 4361 женщин из Ботсваны:
botswana.tsv
О каждой из них мы знаем:
сколько детей она родила (признак ceb)
возраст (age)
длительность получения образования (educ)
религиозная принадлежность (religion)
идеальное, по её мнению, количество детей в семье (idlnchld)
была ли она когда-нибудь замужем (evermarr)
возраст первого замужества (agefm)
длительность получения образования мужем (heduc)
знает ли она о методах контрацепции (knowmeth)
использует ли она методы контрацепции (usemeth)
живёт ли она в городе (urban)
есть ли у неё электричество, радио, телевизор и велосипед (electric, radio, tv, bicycle)
Давайте научимся оценивать количество детей ceb по остальным признакам.
Загрузите данные и внимательно изучите их. Сколько разных значений принимает признак religion?
In [2]:
botswana = pd.read_csv('botswana.tsv', delimiter='\t')
botswana.info()
In [3]:
botswana.columns
botswana.head()
Out[3]:
Out[3]:
In [4]:
print(botswana.religion.value_counts())
Во многих признаках есть пропущенные значения. Сколько объектов из 4361 останется, если выбросить все, содержащие пропуски?
In [5]:
bw_wo_nan = botswana.dropna()
bw_wo_nan.shape
Out[5]:
В разных признаках пропуски возникают по разным причинам и должны обрабатываться по-разному.
Например, в признаке agefm пропуски стоят только там, где evermarr=0, то есть, они соответствуют женщинам, никогда не выходившим замуж. Таким образом, для этого признака NaN соответствует значению "не применимо".
В подобных случаях, когда признак x1 на части объектов в принципе не может принимать никакие значения, рекомендуется поступать так:
создать новый бинарный признак x2={1,0,x1='не применимо',иначе; заменить "не применимо" в x1 на произвольную константу c, которая среди других значений x1 не встречается. Теперь, когда мы построим регрессию на оба признака и получим модель вида y=β0+β1x1+β2x2, на тех объектах, где x1 было измерено, регрессионное уравнение примет вид y=β0+β1x, а там, где x1 было "не применимо", получится y=β0+β1c+β2. Выбор c влияет только на значение и интерпретацию β2, но не β1.
Давайте используем этот метод для обработки пропусков в agefm и heduc.
Создайте признак nevermarr, равный единице там, где в agefm пропуски. Удалите признак evermarr — в сумме с nevermarr он даёт константу, значит, в нашей матрице X будет мультиколлинеарность. Замените NaN в признаке agefm на cagefm=0. У объектов, где nevermarr = 1, замените NaN в признаке heduc на cheduc1=−1 (ноль использовать нельзя, так как он уже встречается у некоторых объектов выборки). Сколько осталось пропущенных значений в признаке heduc?
In [6]:
botswana['nevermarr'] = [1 if botswana.loc[i, 'evermarr'] == 0 else 0 for i in range(botswana.shape[0])]
botswana.head()
Out[6]:
In [7]:
np.unique(botswana.agefm[botswana.agefm.notnull()].values)
Out[7]:
In [8]:
botswana.agefm[botswana.agefm.isnull()] = 0
botswana.head()
Out[8]:
In [9]:
del botswana['evermarr']
In [10]:
botswana.heduc[botswana.heduc.isnull() & botswana.nevermarr.values == 1] = -1
botswana.head()
Out[10]:
In [11]:
botswana.info()
In [12]:
botswana.heduc.isnull().value_counts()
Out[12]:
Избавимся от оставшихся пропусков.
Для признаков idlnchld, heduc и usemeth проведите операцию, аналогичную предыдущей: создайте индикаторы пропусков по этим признакам (idlnchld_noans, heduc_noans, usemeth_noans), замените пропуски на нехарактерные значения (cidlnchld=−1, cheduc2=−2 (значение -1 мы уже использовали), cusemeth=−1).
Остались только пропуски в признаках knowmeth, electric, radio, tv и bicycle. Их очень мало, так что удалите объекты, на которых их значения пропущены.
Какого размера теперь наша матрица данных? Умножьте количество строк на количество всех столбцов (включая отклик ceb).
In [13]:
botswana['idlnchld_noans'] = 0
botswana.loc[botswana.idlnchld.isnull(), 'idlnchld_noans'] = 1
botswana['heduc_noans'] = 0
botswana.loc[botswana.heduc.isnull(), 'heduc_noans'] = 1
botswana['usemeth_noans'] = 0
botswana.loc[botswana.usemeth.isnull(), 'usemeth_noans'] = 1
botswana.head()
Out[13]:
In [16]:
botswana.idlnchld[botswana.idlnchld.isnull()] = -1
botswana.heduc[botswana.heduc.isnull()] = -2
botswana.usemeth[botswana.usemeth.isnull()] = -1
In [17]:
botswana.info()
In [18]:
botswana = botswana.dropna()
In [19]:
botswana.info()
In [20]:
elem_num = botswana.shape[0] * botswana.shape[1]
print('Array size: %d' % elem_num)
Постройте регрессию количества детей ceb на все имеющиеся признаки методом smf.ols, как в разобранном до этого примере. Какой получился коэффициент детерминации R2? Округлите до трёх знаков после десятичной точки.
Если код из примера у вас не воспроизводится:
убедитесь, что вы сделали так:
import statsmodels.formula.api as smf возможно, вам нужно обновить библиотеку patsy; выполните в командной строке
pip install -U patsy
In [21]:
botswana.columns
Out[21]:
In [22]:
formula = 'ceb ~ ' + ' + '.join(botswana.columns[1:])
formula
Out[22]:
In [23]:
reg_m = smf.ols(formula, data=botswana)
fitted_m = reg_m.fit()
print(fitted_m.summary())
In [24]:
print(botswana.religion.value_counts())
Проверьте критерием Бройша-Пагана гомоскедастичность ошибки в построенной модели. Выполняется ли она?
Если ошибка гетероскедастична, перенастройте модель, сделав поправку Уайта типа HC1.
In [25]:
print('Breusch-Pagan test: p=%f' % sms.het_breuschpagan(fitted_m.resid, fitted_m.model.exog)[1])
In [26]:
reg_m2 = smf.ols(formula, data=botswana)
fitted_m2 = reg_m2.fit(cov_type='HC1')
print(fitted_m2.summary())
Удалите из модели незначимые признаки religion, radio и tv. Проверьте гомоскедастичность ошибки, при необходимости сделайте поправку Уайта.
Не произошло ли значимого ухудшения модели после удаления этой группы признаков? Проверьте с помощью критерия Фишера. Чему равен его достигаемый уровень значимости? Округлите до четырёх цифр после десятичной точки.
Если достигаемый уровень значимости получился маленький, верните все удалённые признаки; если он достаточно велик, оставьте модель без религии, тв и радио.
In [32]:
formula2 = 'ceb ~ age + educ + idlnchld + knowmeth + usemeth + agefm + heduc + urban + electric + bicycle \
+ nevermarr + idlnchld_noans + heduc_noans + usemeth_noans'
formula2
Out[32]:
In [33]:
reg_m3 = smf.ols(formula2, data=botswana)
fitted_m3 = reg_m3.fit()
print(fitted_m3.summary())
In [34]:
print('Breusch-Pagan test: p=%f' % sms.het_breuschpagan(fitted_m3.resid, fitted_m3.model.exog)[1])
In [35]:
reg_m4 = smf.ols(formula2, data=botswana)
fitted_m4 = reg_m4.fit(cov_type='HC1')
print(fitted_m4.summary())
In [36]:
print('F=%f, p=%f, k1=%f' % reg_m2.fit().compare_f_test(reg_m4.fit()))
Признак usemeth_noans значим по критерию Стьюдента, то есть, при его удалении модель значимо ухудшится. Но вообще-то отдельно его удалять нельзя: из-за того, что мы перекодировали пропуски в usemeth произвольно выбранным значением cusemeth=−1, удалять usemeth_noans и usemeth можно только вместе.
Удалите из текущей модели usemeth_noans и usemeth. Проверьте критерием Фишера гипотезу о том, что качество модели не ухудшилось. Введите номер первой значащей цифры в достигаемом уровне значимости (например, если вы получили 5.5×10−8, нужно ввести 8).
Если достигаемый уровень значимости получился маленький, верните удалённые признаки; если он достаточно велик, оставьте модель без usemeth и usemeth_noans.
In [37]:
formula3 = 'ceb ~ age + educ + idlnchld + knowmeth + agefm + heduc + urban + electric + bicycle \
+ nevermarr + idlnchld_noans + heduc_noans'
formula3
Out[37]:
In [38]:
reg_m5 = smf.ols(formula3, data=botswana)
fitted_m5 = reg_m5.fit()
print(fitted_m5.summary())
In [44]:
print('F=%f, p=%.40f, k1=%f' % reg_m4.fit().compare_f_test(reg_m5.fit()))
Посмотрите на доверительные интервалы для коэффициентов итоговой модели (не забудьте использовать поправку Уайта, если есть гетероскедастичность ошибки) и выберите правильные выводы.
In [45]:
print(fitted_m4.summary())
In [ ]: